---
title: "MAR Replication"
author: "ASD"
date: "`r Sys.Date()`"
output:
  html_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = F)
knitr::opts_chunk$set(message = F)
knitr::opts_chunk$set(warning = F)

```

```{r}
#load libraries
library(effectsize)
library(equivtest)
library(cregg)
library(broom)
library(multiwayvcov)
library(ggh4x)
library(gridExtra)
library(sandwich)
library(stargazer)
library(tidyverse)
```

```{r}
#load data 
mar_dat <- read_rds("mar_replication_dat.rds") 
### Functions
robust_se <- function(x) {sqrt(diag(vcovHC(x, type = "HC0")))}

cluster_se <- function(x) {sqrt(diag(cluster.vcov(x, mar_dat$respondent)))}
```

## Main Analyses 

### H1 (Fig 2) Does placement of sensitive items (MAR) influence the measurement of those questions? 

```{r}
mar_dat %>% 
  distinct(respondent, .keep_all = T) %>% 
  mutate(placement = if_else(placement == "Two-wave", "Two- \nwave", placement) %>% factor(levels = c("Two- \nwave", "Post", "Pre"))) %>% 
  group_by(placement) %>% 
  dplyr::summarize(avg_mar = mean(mar, na.rm = T)) %>% 
  bind_cols(mar_dat %>% 
    lm(formula = mar ~ placement) %>% 
    cluster_se() %>% 
    tidy()) %>% 
  rename(se = x) %>% 
  mutate(mi = 1.96 * se,
         plus = avg_mar + mi,
         minus = avg_mar - mi) %>%
  ggplot(aes(x = placement, y = avg_mar, color = fct_rev(placement), shape = fct_rev(placement))) +
  geom_point(size = 3) + 
  geom_errorbar(aes(ymin = minus, ymax = plus), width = .2, size= 0.8) + 
  theme_bw() +
  ylim(0.3, 0.42) +
  coord_flip() +
  scale_shape_manual(
    name = "Placement of Sensitive Items",
    values = c(1, 2, 3),
  ) +  
  scale_color_manual(
    name = "Placement of Sensitive Items",
    values = c("grey65", "grey42", "black"),
  ) +
  labs(x = "", y = "Mean Muslim American Resentment") +
    theme(axis.text=element_text(size=13),
        axis.title=element_text(size=13),
        legend.text=element_text(size=9),
        legend.title=element_text(size=11),
        legend.position = "none",
        legend.key.width = unit(1.5,"cm"),
        legend.spacing.x = unit(0.25,"cm"),
        legend.justification = "center") +
  ylim(0,1)
```

#### Corresponding LM (Appendix C.5, Table C.12)

```{r, include = F}
mar_dat %>% 
  distinct(respondent, .keep_all = T) %>% 
  lm(formula = mar ~ placement) %>% 
  stargazer(type = "text", covariate.labels = c("Placement of MAR = Pre", "Placement of MAR = Post"), dep.var.labels = c("Level of MAR"), report=('vc*s*p'), star.cutoffs = c(0.05, 0.01, 0.001))
```

### H2 (Fig 3) Does the placement of the senstive items (MAR) impact green card given?

```{r}
mar_dat %>% 
  mutate(placement = if_else(placement == "Two-wave", "Two- \nwave", placement) %>% factor(levels = c("Two- \nwave", "Post", "Pre"))) %>% 
  group_by(placement) %>% 
  dplyr::summarize(avg_threat = mean(threat_culture, na.rm = T)) %>% 
  bind_cols(mar_dat %>% 
    lm(formula = threat_culture ~ placement) %>% 
    cluster_se() %>% 
    tidy()) %>% 
  rename(se = x) %>% 
  mutate(mi = 1.96 * se,
         plus = avg_threat + mi,
         minus = avg_threat - mi) %>%
  ggplot(aes(x = placement, y = avg_threat, color = fct_rev(placement), shape = fct_rev(placement))) +
  geom_point(size = 3) + 
  geom_errorbar(aes(ymin = minus, ymax = plus), width = .2, size = 0.8) + 
  theme_bw() +
  coord_flip() +
  scale_shape_manual(
    name = "Placement of Sensitive Items",
    values = c(1, 2, 3),
  ) +  
  scale_color_manual(
    name = "Placement of Sensitive Items",
    values = c("grey65", "grey42", "black"),
  ) +
  labs(x = "", y = "Mean Cultural Threat") +
    theme(axis.text=element_text(size=13),
        axis.title=element_text(size=13),
        legend.text=element_text(size=9),
        legend.title=element_text(size=11),
        legend.position = "none",
        legend.key.width = unit(1.5,"cm"),
        legend.spacing.x = unit(0.25,"cm"),
        legend.justification = "center") +
  ylim(0,1)
```

#### Corresponding LM (Appendix C.5, Table C.12)

```{r}
mar_dat %>% 
  lm(formula = threat_culture ~ placement) %>%
  stargazer(type = "text", se = list(cluster_se(.)), covariate.labels = c("Placement of MAR = Pre", "Placement of MAR = Post"), dep.var.labels = c("Cultural Threat"), report=('vc*s*p'), star.cutoffs = c(0.05, 0.01, 0.001))
```

### H3 (Fig 4) Does the placement of sensitive items (MAR) and the substantive treatment impact green card given?

(Set up)
```{r}
#### CREGG Setup to Increase Line Width####
make_feature_headers <- function(x, fmt = "(%s)") {
  feature_levels <- rev(split(x$level, x$feature))
  for (i in seq_along(feature_levels)) {
    feature_levels[[i]] <- levels(x$level)[match(feature_levels[[i]], levels(x$level))]
    feature_levels[[i]] <- c(feature_levels[[i]], sprintf(fmt, names(feature_levels)[i]))
  }
  factor(as.character(x$level), levels = unique(unname(unlist(feature_levels))))
}

plot.cj_mm2 <- 
  function(
    x,
    group = attr(x, "by"),
    feature_headers = TRUE,
    header_fmt = "    (%s)",
    size = 1.0,
    xlab = "Marginal Mean",
    ylab = "",
    legend_title = if (is.null(group)) "Feature" else "Placement of \n Sensitive Items",
    legend_pos = "bottom",
    xlim = NULL,
    vline = 0,
    vline_color = "gray",
    theme = ggplot2::theme_bw(),
    ...
  ) {
    
    # optionally, add gaps between features
    if (isTRUE(feature_headers)) {
      x$level <- make_feature_headers(x, fmt = header_fmt)
      to_merge <- data.frame(feature = unique(x$feature), level = sprintf(header_fmt, unique(x$feature)))
      if ("BY" %in% names(x)) {
        to_merge <- do.call("rbind", lapply(unique(x[["BY"]]), function(lev) {
          to_merge[["BY"]] <- lev
          to_merge
        }))
      } else if (!is.null(group)) {
        to_merge <- do.call("rbind", lapply(unique(x[[group]]), function(lev) {
          to_merge[[group]] <- lev
          to_merge
        }))
      }
      x <- merge(x, to_merge, all = TRUE)
    }
    
    if (is.null(group)) {
      p <- ggplot2::ggplot(data = x, ggplot2::aes_string(x = "estimate", y = "level", colour = "feature"))
    } else {
      if (is.null(x[[group]])) {
        stop(sprintf("`group` variable '%s' not found", group))
      }
      p <- ggplot2::ggplot(data = x, ggplot2::aes_string(x = "estimate", y = "level", colour = group, group = group))
    }
    
    if (is.null(xlim)) {
      xmin <- min(x$lower, na.rm = TRUE)
      xmin <- if (xmin < 0) 1.04*xmin else .96*xmin
      xmax <- max(x$upper, na.rm = TRUE)
      xmax <- if (xmax > 0) 1.04*xmax else .96*xmax
      xlim <- c(xmin, xmax)
    }
    
    if (!is.null(vline)) {
      p <- p + ggplot2::geom_vline(xintercept = vline, colour = vline_color)
    }
    
    p <- p + ggplot2::geom_point(position = ggstance::position_dodgev(height = 0.75), size = size, na.rm = TRUE) +
      ggplot2::geom_errorbarh(ggplot2::aes_string(xmin = "lower", xmax = "upper"),  
                              size = 0.6, height = 0, na.rm = TRUE,
                              position = ggstance::position_dodgev(height = 0.75))
    if (is.null(group)) {
      p <- p + ggplot2::scale_colour_discrete(guide = ggplot2::guide_legend(title = legend_title))
    } else {
      p <- p + ggplot2::scale_colour_discrete(breaks = levels(x[[group]]),
                                              labels = levels(x[[group]]),
                                              guide = ggplot2::guide_legend(title = legend_title))
    }
    p <- p +
      ggplot2::scale_x_continuous(limits = xlim, oob = scales::rescale_none) +
      ggplot2::xlab(xlab) + 
      ggplot2::ylab(ylab) + 
      theme + ggplot2::theme(
        legend.position = legend_pos,
        panel.grid.major = ggplot2::element_blank(),
        panel.grid.minor = ggplot2::element_blank()
      ) + 
      ggplot2::guides(colour = ggplot2::guide_legend(title = legend_title))
    return(p)
  }

```

```{r}
mar_dat %>%
  mutate(
    placement = placement %>% factor(levels = c("Two-wave", "Post", "Pre")),
    race = race %>% as.factor() %>% factor(levels = c(
      "White", "Black", "Middle Eastern", "South Asian"
    ))
  ) %>%
  rename(Religion = religion,
         Race = race)  %>%
  cj(
    formula = greencard ~ Religion + Race,
    id = ~ respondent,
    by = ~ placement,
    estimate = c("mm")
  ) %>%
  mutate(level2 = level %>% factor(
    levels = c(
      "White",
      "Black",
      "Middle Eastern",
      "South Asian",
      "    (Race)",
      "Christian",
      "Jewish",
      "Muslim",
      "    (Religion)"
    )
  )) %>%
  ggplot(aes(
    x = estimate,
    y = level2,
    shape = placement,
    color = placement
  )) +
  geom_point(position = position_dodge(0.7), size = 2) +
  geom_line(position = position_dodge(0.7), size = 1) +
  geom_errorbar(aes(xmin = lower, xmax = upper),
                width = 0.65,
                size = 0.65,
                position = position_dodge(0.7), 
                size = 1) +
  geom_vline(xintercept = 0.5, color = "grey50")+
  scale_color_manual(
    "Placement of Sensitive Items",
    values = rev(c("grey65", "grey42", "black")),
    drop = FALSE
  ) +
  scale_shape_manual(
    name = "Placement of Sensitive Items",
    values = c(3, 2, 1)
  ) +
  labs(y = "", x = "Marginal Mean Green Card Given") +
  scale_y_discrete(drop = FALSE) +
  xlim(0.35, 0.65) +
  theme_bw() +
  theme(
    axis.text=element_text(size=13),
        axis.title=element_text(size=13),
        legend.text=element_text(size=9),
        legend.title=element_text(size=11),
    legend.position = "bottom",
    legend.key.width = unit(1,"cm"),
    legend.spacing.x = unit(0.3,"cm"),
    legend.box.background = element_rect(colour = "black"),
    legend.justification = "center",
    axis.text.y = element_text(hjust = (0)),
    axis.text.x = element_text(face = "plain"),
    axis.ticks = element_blank(),
    text = element_text(
      size = 12,
      face = c(
        #race
        "plain",
        #White
        "plain",
        #Black
        "plain",
        #ME
        "plain",
        #SA
        "bold",
        #race
        #religion
        "plain",
        #christian
        "plain",
        #jewish
        "plain",
        #muslim
        "bold"
        #religion      
        )
    )  )  +
  guides(shape = guide_legend(title.position="top", 
                              title.hjust = 0.5,
                              reverse = T), 
        color = guide_legend(title.position="top", 
                              title.hjust = 0.5,
                              reverse = T)
        ) +
  xlim(0,1)
```

#### H3 LMs (Appendix C.5, Table C.13)
```{r}
mar_dat %>%
  lm(formula = greencard ~ religion * placement) %>%
  stargazer(
    type = "text",
    se = list(cluster_se(.)),
    dep.var.labels = c("Green Card to Immigrant"),
    star.cutoffs = c(0.05, 0.01, 0.001)
    )
```

#### H3 LMs (Appendix C.5, Table C.14)
```{r}
mar_dat %>%
  lm(formula = greencard ~ race * placement) %>%
  stargazer(
    type = "text",
    se = list(cluster_se(.)),
    dep.var.labels = c("Green Card to Immigrant"),
    star.cutoffs = c(0.05, 0.01, 0.001) 
    )
```

#### H3 Difference pre-/post-treatment and two-wave (Appendix C.5, Table C.15)
```{r}
#Rows 1-14
mar_dat %>%
  rename(
    Religion = religion,
    Race = race
  ) %>%
  cj(
    formula = greencard ~ Religion  + Race,
    id = ~ respondent,
    by = ~ placement,
    estimate = c("mm_diff")
  ) %>%
  select(BY, feature, level, estimate, std.error, p) %>%
  rename(Attribute = feature, 
         Level = level, 
         SE = std.error,
         `Marginal Mean` = estimate) %>% 
  mutate(p.adj = p.adjust(p, "fdr")) %>% 
  knitr::kable(digits = 3)

#Rows 15-21
mar_dat %>%
  mutate(placement = placement %>% factor(levels = c("Pre", "Post", "Two-wave"))) %>% 
  rename(
    Religion = religion,
    Race = race
  ) %>%
  cj(
    formula = greencard ~ Religion  + Race,
    id = ~ respondent,
    by = ~ placement,
    estimate = c("mm_diff")
  ) %>%
  slice(1:7) %>% 
  select(BY, feature, level, estimate, std.error, p) %>%
  rename(Attribute = feature, 
         Level = level, 
         SE = std.error,
         `Marginal Mean` = estimate) %>% 
    mutate(p.adj = p.adjust(p, "fdr")) %>% 
   knitr::kable(digits = 3)

# Custom function to calculate p-values
calculate_p_value <- function(estimate, se) {
  t_statistic <- estimate / se
  p_value <- 2 * (1 - pt(abs(t_statistic), 5925))
  return(p_value)
}

#Rows 22-42
mar_dat %>%
  mutate(placement = placement %>%  factor(levels = c("Pre", "Post", "Two-wave"))) %>% 
  rename(
    Religion = religion,
    Race = race
  ) %>%
  cj(
    formula = greencard ~ Religion  + Race,
    id = ~ respondent,
    by = ~ placement,
    estimate = c("mm")
  ) %>% 
  mutate(est = 0.5 - estimate, # p value has to be relative to 0.5 as the baseline 
         pv = map2(est, std.error, calculate_p_value),
         pv = as.numeric(pv)
    ) %>% 
  select(BY, estimate, std.error, pv) %>% 
  knitr::kable(digits = 3)
```


##### H3 Full MM Graph (Appendix C.5, Figure C.2)

```{r}
mar_dat %>%
  mutate(
    education = factor(
      education,
      levels = c("Elementary school", "High school", "College", "Master's Degree")
    ),
    english = factor(english, levels = c("Intermediate", "Advanced", "Fluent")),
    race = race %>% as.factor() %>% factor(levels = c(
      "White", "Black", "Middle Eastern", "South Asian"
    ))
  ) %>%
  rename(
    Education = education,
    Gender = gender,
    `English Fluency` = english,
    Religion = religion,
    Race = race
  ) %>%
  cj(
    formula = greencard ~ Education + Gender + Religion + `English Fluency` + Race,
    id = ~ respondent,
    by = ~ placement,
    estimate = c("mm")
  ) %>%
  filter(!is.na(placement)) %>% 
  plot.cj_mm2(size = 2,
              vline = 0.5, group = "placement") +
  aes(shape = placement) +
  scale_color_manual(
    name = "Placement of \n Sensitive Items",
    values = rev(c("grey65", "grey42", "black")),
    na.translate = F,
  ) +
scale_shape_manual(
    name = "Placement of \n Sensitive Items",
    values = c(3,2,1),
    na.translate = F,
  ) +
  theme(
    axis.text=element_text(size=11),
        axis.title=element_text(size=12),
        legend.text=element_text(size=9),
        legend.title=element_text(size=11),
        legend.position = c(0.85, 0.20),
                       legend.spacing.y = unit(1, "mm"),
                       panel.border = element_rect(colour = "black", fill=NA),
                       legend.background = element_blank(),
                       legend.box.background = element_rect(colour = "black"),
    axis.text.y = element_text(hjust = (0)),
    axis.text.x = element_text(face = "plain"),
    axis.ticks = element_blank(),
    text = element_text(
      size = 16,
      face = c(
        #race
        "plain",
        #White
        "plain",
        #Black
        "plain",
        #ME
        "plain",
        #SA
        "bold",
        #race
        #religion
        "plain",
        #christian
        "plain",
        #jewish
        "plain",
        #muslimx
        "bold",
        #religion
        #english
        "plain",
        #intermediate
        "plain",
        #advanced
        "plain",
        #fluent
        "bold",
        #english
        #gender
        "plain",
        #Male
        "plain",
        #Female
        "bold",
        #Gender
        #education
        "plain",
        #elementary
        "plain",
        #high school
        "plain",
        #college
        "plain",
        #masters
        "bold" #education
      )
    ),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    axis.line = element_line(colour = "black")
  ) + guides(shape = guide_legend(reverse = T)) + 
  guides(color = guide_legend(reverse = T)) 
```




#### H3 Alternative DV (Cultural Threat) x Religion (Appendix C.5, Table C.16)
```{r}
mar_dat %>%
  lm(formula = threat_culture ~ religion * placement) %>%
  stargazer(
    type = "text",
    se = list(cluster_se(.)),
    dep.var.labels = c("Green Card to Immigrant"),
    star.cutoffs = c(0.05, 0.01, 0.001))

```

#### H3 Alternative DV (Cultural Threat) x Race (Appendix C.5, Table C.17)

```{r}
mar_dat %>%
  lm(formula = threat_culture ~ race * placement) %>%
  stargazer(
    type = "text",
    se = list(cluster_se(.)),
    dep.var.labels = c("Green Card to Immigrant"),
    star.cutoffs = c(0.05, 0.01, 0.001))
```


### H4 (Fig 5) Does the placement of sensitive items (MAR), the measurement of the sensitive items (MAR), and the substantive treatment impact green card given?


```{r}
mar_dat %>%
  mutate(
    mar_hl = if_else(mar < 0.5, "Low MAR", "High MAR") %>% as.factor() %>% factor(levels = c("Low MAR", "High MAR")),
    placement = placement %>% factor(levels = c("Two-wave", "Post", "Pre")),
    race = race %>% as.factor() %>% factor(levels = c(
      "White", "Black", "Middle Eastern", "South Asian"
    ))
  ) %>%
  rename(Religion = religion,
         Race = race)  %>%
  cj(
    formula = greencard ~ Religion + Race,
    id = ~ respondent,
    by = ~ placement * mar_hl,
    estimate = c("mm")
  ) %>%
  mutate(level2 = level %>% factor(
    levels = c(
      "White",
      "Black",
      "Middle Eastern",
      "South Asian",
      "    (Race)",
      "Christian",
      "Jewish",
      "Muslim",
      "    (Religion)"
    )
  )) %>%
  ggplot(aes(
    x = estimate,
    y = level2,
    shape = placement,
    color = placement
  )) +
  geom_point(position = position_dodge(0.7), size = 2) +
  geom_line(position = position_dodge(0.7), size = 1) +
  geom_errorbar(aes(xmin = lower, xmax = upper),
                width = .65,
                size = 0.65,
                position = position_dodge(0.7)) +
  geom_vline(xintercept = 0.5, color = "grey50", alpha = 0.5)+
  scale_color_manual(
    "Placement of Sensitive Items",
    values = rev(c("grey65", "grey42", "black")),
  ) +
  scale_shape_manual(
    name = "Placement of Sensitive Items",
    values = rev(c(1, 2, 3))
  ) +
  labs(y = "", x = "Marginal Mean Green Card Given") +
  scale_y_discrete(drop = FALSE) +
  facet_wrap(~mar_hl) +
  xlim(0.17, 0.75)+
  theme_bw() +
  theme(
    axis.text=element_text(size=13),
    axis.title=element_text(size=13),
    legend.text=element_text(size=9),
    legend.title=element_text(size=11),
    legend.position = "bottom",
    legend.key.width = unit(1,"cm"),
    legend.spacing.x = unit(0.3,"cm"),
    legend.justification = "center",
    legend.box.background = element_rect(colour = "black"),
    axis.text.y = element_text(hjust = (0)),
    axis.text.x = element_text(face = "plain"),
    axis.ticks = element_blank(),
    text = element_text(
      size = 12,
      face = c(
        #race
        "plain",
        #White
        "plain",
        #Black
        "plain",
        #ME
        "plain",
        #SA
        "bold",
        #race
        #religion
        "plain",
        #christian
        "plain",
        #jewish
        "plain",
        #muslim
        "bold"
        #religion      
        )
    ), 
    panel.spacing.x = unit(1, "lines")) +
  guides(shape = guide_legend(title.position="top", 
                              title.hjust = 0.5,
                              reverse = T), 
        color = guide_legend(title.position="top", 
                              title.hjust = 0.5,
                              reverse = T)
        ) +
  xlim(0,1)
```


#### Corresponding LMs
##### Religion (Appendix C.5, Table C.18)

```{r}
mar_dat %>% 
  lm(formula = greencard ~ placement*mar*religion) %>%
  stargazer(se = list(cluster_se(.)), 
                type = "text",
            covariate.labels = c("Pre", "Post", "Muslim American Resentment (MAR)", "Jewish", "Muslim", "Pre x MAR", "Post x MAR", "Pre x Jewish", "Post x Jewish", "Pre x Muslim", "Post x Muslism", "MAR x Jewish", "MAR x Muslim", "Pre x MAR x Jewish", "Post x MAR x Jewish", "Pre x MAR x Muslim", "Post x MAR x Muslim", "Constant"), dep.var.labels = c("Green Card to Immigrant"),
            star.cutoffs = c(0.05, 0.01, 0.001))
```


##### Race (Appendix C.5, Table C.19)
```{r}
mar_dat %>% 
  lm(formula = greencard ~ placement*mar*race) %>%
  stargazer(se = list(cluster_se(.)), 
            type = "text",
            covariate.labels = c("Pre", "Post", "Muslim American Resentment (MAR)", "Black", "Middle Eastern", "South Asian", "Pre x MAR", "Post x MAR", "Pre x Black", "Post x Black", "Pre x Middle Eastern", "Post x Middle Eastern", "Pre x South Asian", "Post x South Asian", "MAR x Black", "MAR x Middle Eastern", "MAR x South Asian", "Pre x MAR x Black", "Post x MAR x Black", "Pre x MAR x Middle Eastern", "Post x MAR x Middle Eastern", "Pre x MAR x South Asian", "Post x MAR x South Asian", "Constant"), dep.var.labels = c("Green Card to Immigrant"), star.cutoffs = c(0.05, 0.01, 0.001))
```

#### H4: Difference in means for Low vs High on MAR
##### Low Mar (Appendix C.5, Table C.20)

```{r}
#Rows 1-14
mar_dat %>% 
  filter(mar < 0.5) %>% 
  mutate(
    placement = placement %>% factor(levels = c("Pre", "Post", "Two-wave"))
  ) %>%
  filter(!is.na(race)) %>%
  rename(
    Religion = religion,
    Race = race
  ) %>%
  cj(
    formula = greencard ~ Religion + Race,
    id = ~ respondent,
    by = ~ placement,
    estimate = c("mm_diff")
  ) %>%
  select(BY, feature, level, estimate, std.error, p) %>%
  rename(Attribute = feature, 
         Level = level, 
         SE = std.error,
         `Marginal Mean` = estimate) %>% 
   mutate(p.adj = p.adjust(p, "fdr")) %>% 
   knitr::kable(digits = 3)

#Rows 14-21
mar_dat %>% 
  filter(mar < 0.5 & placement != "Pre") %>% 
  mutate(
    placement = placement %>% factor(levels = c("Post", "Two-wave"))
  ) %>%
  filter(!is.na(race)) %>%
  rename(
    Religion = religion,
    Race = race
  ) %>%
  cj(
    formula = greencard ~ Religion + Race,
    id = ~ respondent,
    by = ~ placement,
    estimate = c("mm_diff")
  ) %>%
  select(BY, feature, level, estimate, std.error, p) %>%
  rename(Attribute = feature, 
         Level = level, 
         SE = std.error,
         `Marginal Mean` = estimate)  %>%
  mutate(p.adj = p.adjust(p, "fdr")) %>% 
  knitr::kable(digits = 3)

# Custom function to calculate p-values
calculate_p_value <- function(estimate, se) {
  t_statistic <- estimate / se
  p_value <- 2 * (1 - pt(abs(t_statistic), 4307))
  return(p_value)
}

#Rows 22-42
mar_dat %>%
  rename(Religion = religion, 
         Race = race) %>% 
  mutate(placement = placement %>% factor(levels = c("Pre", "Post", "Two-wave"))) %>% 
  filter(mar < 0.5) %>% 
  cj(
    formula = greencard ~ Religion  + Race,
    id = ~ respondent,
    by = ~ placement,
    estimate = c("mm")
  ) %>% 
  mutate(est = estimate - 0.5, # p value has to be relative to 0.5 as the baseline 
         pv = map2(est, std.error, calculate_p_value), 
         pv = as.numeric(pv),
         p.adj = p.adjust(pv, "fdr")
    ) %>% 
  select(BY, feature, level, estimate, std.error, pv, p.adj) %>% 
  knitr::kable(digit = 3)
```

##### High Mar (Appendix C.5, Table C.20 Cont.)

```{r}
#Rows 43-56
mar_dat %>% 
  filter(mar >= 0.5) %>% 
  mutate(
    placement = placement %>% factor(levels = c("Pre", "Post", "Two-wave"))
  ) %>%
  rename(
    Religion = religion,
    Race = race
  ) %>%
  cj(
    formula = greencard ~ Religion + Race,
    id = ~ respondent,
    by = ~ placement,
    estimate = c("mm_diff")
  ) %>%
  select(BY, feature, level, estimate, std.error, p) %>%
  rename(Attribute = feature, 
         Level = level, 
         SE = std.error,
         `Marginal Mean` = estimate) %>% 
   mutate(p.adj = p.adjust(p, "fdr")) %>% 
   knitr::kable(digits = 3)

#Rows 57-63
mar_dat %>% 
  filter(mar >= 0.5 & placement != "Pre") %>%
  mutate(placement = placement %>% factor(levels = c("Post", "Two-wave"))) %>%
  filter(!is.na(race)) %>%
  rename(
    Religion = religion,
    Race = race
  ) %>%
  cj(
    formula = greencard ~ Religion + Race,
    id = ~ respondent,
    by = ~ placement,
    estimate = c("mm_diff")
  ) %>%
  select(BY, feature, level, estimate, std.error, p) %>%
  rename(Attribute = feature, 
         Level = level, 
         SE = std.error,
         `Marginal Mean` = estimate) %>%
  mutate(p.adj = p.adjust(p, "fdr")) %>% 
   knitr::kable(digits = 3)

# Custom function to calculate p-values
calculate_p_value <- function(estimate, se) {
  t_statistic <- estimate / se
  p_value <- 2 * (1 - pt(abs(t_statistic), 1606))
  return(p_value)
}

#Rows 63-84
mar_dat %>%
  mutate(placement = placement %>% factor(levels = c("Pre", "Post", "Two-wave"))) %>% 
  filter(mar >= 0.5) %>% 
  cj(
    formula = greencard ~ religion  + race,
    id = ~ respondent,
    by = ~ placement,
    estimate = c("mm")
  ) %>% 
  mutate(est = estimate - 0.5, # p value has to be relative to 0.5 as the baseline 
         pv = map2(est, std.error, calculate_p_value), 
         pv = as.numeric(pv), 
         p.adj = p.adjust(pv, "fdr")
    ) %>% 
  select(BY, feature, level, estimate, std.error, pv, p.adj) %>% 
  knitr::kable(digits = 3)

```

#### Full H4 (Appendic C.5, Figure C.3)
```{r}
mar_dat %>%
  mutate(
    mar_hl = if_else(mar < 0.5, "Low MAR", "High MAR") %>% as.factor() %>% factor(levels = c("Low MAR", "High MAR")),
    education = factor(
      education,
      levels = c("Elementary school", "High school", "College", "Master's Degree")
    ),
    english = factor(english, levels = c("Intermediate", "Advanced", "Fluent")),
    placement = placement %>% factor(levels = c("Two-wave", "Post", "Pre")),
    race = race %>% as.factor() %>% factor(levels = c(
      "White", "Black", "Middle Eastern", "South Asian"
    ))) %>%
  rename(
    Education = education,
    Gender = gender,
    `English Fluency` = english,
    Religion = religion,
    Race = race
  ) %>%
  cj(
    formula = greencard ~ Education + Gender + Religion + `English Fluency` + Race,
    id = ~ respondent,
    by = ~ placement * mar_hl,
    estimate = c("mm")
  ) %>%
  mutate(level2 = level %>% factor(
    levels = c(
      "White",
      "Black",
      "Middle Eastern",
      "South Asian",
      "    (Race)",
      "Intermediate",
      "Advanced",
      "Fluent",
      "    (English Fluencey)",
      "Christian",
      "Jewish",
      "Muslim",
      "    (Religion)",
      "Male",
      "Female",
      "    (Gender)",
      "Elementary school",
      "High school",
      "College",
      "Master's Degree",
      "    (Education)"
    )
  )) %>%
  ggplot(aes(
    x = estimate,
    y = level2,
    shape = placement,
    color = placement
  )) +
  geom_point(position = position_dodge(0.7), size = 2) +
  geom_line(position = position_dodge(0.7), size = 1) +
  geom_errorbar(aes(xmin = lower, xmax = upper),
                width = .65,
                size = 0.65,
                position = position_dodge(0.7)) +
  geom_vline(xintercept = 0.5, color = "grey50", alpha = 0.5)+
  scale_color_manual(
    "Placement of Sensitive Items",
    values = rev(c("grey65", "grey42", "black")),
  ) +
  scale_shape_manual(
    name = "Placement of Sensitive Items",
    values = rev(c(1, 2, 3))
  ) +
  labs(y = "", x = "Marginal Mean Green Card Given") +
  scale_y_discrete(drop = FALSE) +
  facet_wrap(~mar_hl) +
  theme_bw() +
  theme(
    axis.text=element_text(size=13),
    axis.title=element_text(size=13),
    legend.text=element_text(size=9),
    legend.title=element_text(size=11),
    legend.position = "bottom",
    legend.key.width = unit(1,"cm"),
    legend.spacing.x = unit(0.3,"cm"),
    legend.justification = "center",
    legend.box.background = element_rect(colour = "black"),
    axis.text.y = element_text(hjust = (0)),
    axis.text.x = element_text(face = "plain"),
    axis.ticks = element_blank(),
    text = element_text(
      size = 12,
      face = c(
        #race
        "plain",
        #White
        "plain",
        #Black
        "plain",
        #ME
        "plain",
        #SA
        "bold",
        #race
        #religion
        "plain",
        #christian
        "plain",
        #jewish
        "plain",
        #muslimx
        "bold",
        #religion
        #english
        "plain",
        #intermediate
        "plain",
        #advanced
        "plain",
        #fluent
        "bold",
        #english
        #gender
        "plain",
        #Male
        "plain",
        #Female
        "bold",
        #Gender
        #education
        "plain",
        #elementary
        "plain",
        #high school
        "plain",
        #college
        "plain",
        #masters
        "bold" #education
      )
    ) )+
  guides(shape = guide_legend(title.position="top", 
                              title.hjust = 0.5,
                              reverse = T), 
        color = guide_legend(title.position="top", 
                              title.hjust = 0.5,
                              reverse = T)
        )
    
#ggsave("appendix_mar4.png", width = 250, height = 125, units = "mm", dpi=1200)
```

#### Alternative DV: Cultural Threat 
##### Religion (Appendic C.5, Table C.21)
```{r}
mar_dat %>% 
  lm(formula = threat_culture ~ placement*mar*religion) %>%
  stargazer(se = list(cluster_se(.)), 
                type = "text",
            covariate.labels = c("Pre", "Post", "Muslim American Resentment (MAR)", "Jewish", "Muslim", "Pre x MAR", "Post x MAR", "Pre x Jewish", "Post x Jewish", "Pre x Muslim", "Post x Muslism", "MAR x Jewish", "MAR x Muslim", "Pre x MAR x Jewish", "Post x MAR x Jewish", "Pre x MAR x Muslim", "Post x MAR x Muslim", "Constant"), dep.var.labels = c("Cultural Threat"), star.cutoffs = c(0.05, 0.01, 0.001))
```

##### Race (Appendix C.5, Table C.22)
```{r}
mar_dat %>% 
  lm(formula = threat_culture ~ placement*mar*race) %>%
  stargazer(se = list(cluster_se(.)), 
            type = "text",
            covariate.labels = c("Pre", "Post", "Muslim American Resentment (MAR)", "Black", "Middle Eastern", "South Asian", "Pre x MAR", "Post x MAR", "Pre x Black", "Post x Black", "Pre x Middle Eastern", "Post x Middle Eastern", "Pre x South Asian", "Post x South Asian", "MAR x Black", "MAR x Middle Eastern", "MAR x South Asian", "Pre x MAR x Black", "Post x MAR x Black", "Pre x MAR x Middle Eastern", "Post x MAR x Middle Eastern", "Pre x MAR x South Asian", "Post x MAR x South Asian", "Constant"), dep.var.labels = c("Cultural Threat"), star.cutoffs = c(0.05, 0.01, 0.001)
           )
```




## Equivalence Test (Appendix C.5, Table C.23)

```{r}
# H1: Placement of sensitive item on measurement of that sensitive item
data1 <- mar_dat$mar[mar_dat$placement == "Pre"]
data2 <- mar_dat$mar[mar_dat$placement == "Post"]
data3 <- mar_dat$mar[mar_dat$placement == "Two-wave"]


equiv.t.test(data1, data2, eps_std = 0.11) %>% summary()
equiv.t.test(data1, data3, eps_std = 0.11) %>% summary()
equiv.t.test(data2, data3, eps_std = 0.17) %>% summary()


# H2: Placement of sensitive item on measure of experimental DV 
data1 <- mar_dat$greencard[mar_dat$placement == "Pre"]
data2 <- mar_dat$greencard[mar_dat$placement == "Post"]
data3 <- mar_dat$greencard[mar_dat$placement == "Two-wave"]

equiv.t.test(data1, data2, eps_std = 0.01) %>% summary()
equiv.t.test(data1, data3, eps_std = 0.01) %>% summary()
equiv.t.test(data2, data3, eps_std = 0.01) %>% summary()

# H3
m2 <- mar_dat %>% 
  unite("race", placement, race, remove = F) %>% 
  unite("religion", placement, religion)

equiv.t.test(m2$greencard[m2$religion == "Pre_Christian"], m2$greencard[m2$religion == "Post_Christian"], eps_std = 0.13) %>% summary()
equiv.t.test(m2$greencard[m2$religion == "Pre_Jewish"], m2$greencard[m2$religion ==    "Post_Jewish"], eps_std = 0.14) %>% summary() 
equiv.t.test(m2$greencard[m2$religion == "Pre_Muslim"], m2$greencard[m2$religion ==    "Post_Muslim"], eps_std = 0.19) %>% summary() 

equiv.t.test(m2$greencard[m2$race == "Pre_White"], m2$greencard[m2$race == "Post_White"], eps_std = 0.19) %>% summary() 
equiv.t.test(m2$greencard[m2$race == "Pre_Black"], m2$greencard[m2$race == "Post_Black"], eps_std = 0.27) %>% summary() 
equiv.t.test(m2$greencard[m2$race == "Pre_Middle Eastern"], m2$greencard[m2$race == "Post_Middle Eastern"], eps_std = 0.13) %>% summary() 
equiv.t.test(m2$greencard[m2$race == "Pre_South Asian"], m2$greencard[m2$race == "Post_South Asian"], eps_std = 0.2) %>% summary() 

equiv.t.test(m2$greencard[m2$religion == "Pre_Christian"], m2$greencard[m2$religion == "Two-wave_Christian"] , eps_std = 0.11) %>% summary()
equiv.t.test(m2$greencard[m2$religion == "Pre_Jewish"], m2$greencard[m2$religion == "Two-wave_Jewish"], eps_std = 0.09) %>% summary() 
equiv.t.test(m2$greencard[m2$religion == "Pre_Muslim"], m2$greencard[m2$religion == "Two-wave_Muslim"], eps_std = 0.01) %>% summary() 

equiv.t.test(m2$greencard[m2$race == "Pre_White"], m2$greencard[m2$race == "Two-wave_White"], eps_std = 0.14) %>% summary() 
equiv.t.test(m2$greencard[m2$race == "Pre_Black"], m2$greencard[m2$race == "Two-wave_Black"], eps_std = 0.23) %>% summary() 
equiv.t.test(m2$greencard[m2$race == "Pre_Middle Eastern"], m2$greencard[m2$race == "Two-wave_Middle Eastern"], eps_std = 0.12) %>% summary() 
equiv.t.test(m2$greencard[m2$race == "Pre_South Asian"], m2$greencard[m2$race == "Two-wave_South Asian"], eps_std = 0.15) %>% summary() 

equiv.t.test(m2$greencard[m2$religion == "Post_Christian"], m2$greencard[m2$religion == "Two-wave_Christian"], eps_std = 0.15) %>% summary() 
equiv.t.test(m2$greencard[m2$religion == "Post_Jewish"], m2$greencard[m2$religion == "Two-wave_Jewish"], eps_std = 0.14) %>% summary() 
equiv.t.test(m2$greencard[m2$religion == "Post_Muslim"], m2$greencard[m2$religion == "Two-wave_Muslim"]  , eps_std = 0.19) %>% summary() 
      
equiv.t.test(m2$greencard[m2$race == "Post_White"], m2$greencard[m2$race == "Two-wave_White"], eps_std = 0.17) %>% summary() 
equiv.t.test(m2$greencard[m2$race == "Post_Black"], m2$greencard[m2$race == "Two-wave_Black"], eps_std = 0.2) %>% summary() 
equiv.t.test(m2$greencard[m2$race == "Post_Middle Eastern"], m2$greencard[m2$race == "Two-wave_Middle Eastern"], eps_std = 0.1) %>% summary() 
equiv.t.test(m2$greencard[m2$race == "Post_South Asian"], m2$greencard[m2$race == "Two-wave_South Asian"], eps_std = 0.16) %>% summary() 

# H4
#Low
m3 <- m2 %>% 
  mutate(mar_hl = if_else(mar < 0.5, 0, 1)) %>% 
  filter(mar_hl == 0)

equiv.t.test(m3$greencard[m3$religion == "Pre_Christian"], m3$greencard[m3$religion == "Post_Christian"], eps_std = 0.14) %>% summary() 
equiv.t.test(m3$greencard[m3$religion == "Pre_Jewish"], m3$greencard[m3$religion == "Post_Jewish"], eps_std = 0.12) %>% summary() 
equiv.t.test(m3$greencard[m3$religion == "Pre_Muslim"], m3$greencard[m3$religion == "Post_Muslim"], eps_std = 0.17) %>% summary() 

equiv.t.test(m3$greencard[m3$race == "Pre_White"], m3$greencard[m3$race == "Post_White"], eps_std = 0.18) %>% summary() 
equiv.t.test(m3$greencard[m3$race == "Pre_Black"], m3$greencard[m3$race == "Post_Black"], eps_std = 0.33) %>% summary() 
equiv.t.test(m3$greencard[m3$race == "Pre_Middle Eastern"], m3$greencard[m3$race == "Post_Middle Eastern"], eps_std = 0.14) %>% summary() 
equiv.t.test(m3$greencard[m3$race == "Pre_South Asian"], m3$greencard[m3$race == "Post_South Asian"], eps_std = 0.21) %>% summary() 

equiv.t.test(m3$greencard[m3$religion == "Pre_Christian"], m3$greencard[m3$religion == "Two-wave_Christian"], eps_std = 0.19) %>% summary() 
equiv.t.test(m3$greencard[m3$religion == "Pre_Jewish"], m3$greencard[m3$religion == "Two-wave_Jewish"], eps_std = 0.16) %>% summary() 
equiv.t.test(m3$greencard[m3$religion == "Pre_Muslim"], m3$greencard[m3$religion == "Two-wave_Muslim"], eps_std = 0.13) %>% summary() 

equiv.t.test(m3$greencard[m3$race == "Pre_White"], m3$greencard[m3$race == "Two-wave_White"], eps_std = 0.12) %>% summary() 
equiv.t.test(m3$greencard[m3$race == "Pre_Black"], m3$greencard[m3$race == "Two-wave_Black"], eps_std = 0.34) %>% summary() 
equiv.t.test(m3$greencard[m3$race == "Pre_Middle Eastern"], m3$greencard[m3$race == "Two-wave_Middle Eastern"], eps_std = 0.16) %>% summary() 
equiv.t.test(m3$greencard[m3$race == "Pre_South Asian"], m3$greencard[m3$race == "Two-wave_South Asian"], eps_std = 0.04) %>% summary() 

equiv.t.test(m3$greencard[m3$religion == "Post_Christian"], m3$greencard[m3$religion == "Two-wave_Christian"], eps_std = 0.16) %>% summary() 
equiv.t.test(m3$greencard[m3$religion == "Post_Jewish"], m3$greencard[m3$religion == "Two-wave_Jewish"], eps_std = 0.19) %>% summary() 
equiv.t.test(m3$greencard[m3$religion == "Post_Muslim"], m3$greencard[m3$religion == "Two-wave_Muslim"]  , eps_std = 0.15) %>% summary() 

equiv.t.test(m3$greencard[m3$race == "Post_White"], m3$greencard[m3$race == "Two-wave_White"], eps_std = 0.17) %>% summary() 
equiv.t.test(m3$greencard[m3$race == "Post_Black"], m3$greencard[m3$race == "Two-wave_Black"], eps_std = 0.15) %>% summary() 
equiv.t.test(m3$greencard[m3$race == "Post_Middle Eastern"], m3$greencard[m3$race == "Two-wave_Middle Eastern"]  , eps_std = 0.13) %>% summary() 
equiv.t.test(m3$greencard[m3$race == "Post_South Asian"], m3$greencard[m3$race == "Two-wave_South Asian"]  , eps_std = .22) %>% summary() 


#High
m3 <- m2 %>% 
  mutate(mar_hl = if_else(mar < 0.5, 0, 1)) %>% 
  filter(mar_hl == 1)

equiv.t.test(m3$greencard[m3$religion == "Pre_Christian"], m3$greencard[m3$religion == "Post_Christian"], eps_std = 0.21) %>% summary() 
equiv.t.test(m3$greencard[m3$religion == "Pre_Jewish"], m3$greencard[m3$religion == "Post_Jewish"], eps_std = 0.35) %>% summary() 
equiv.t.test(m3$greencard[m3$religion == "Pre_Muslim"], m3$greencard[m3$religion == "Post_Muslim"], eps_std = 0.42) %>% summary() 

equiv.t.test(m3$greencard[m3$race == "Pre_White"], m3$greencard[m3$race == "Post_White"], eps_std = 0.41) %>% summary() 
equiv.t.test(m3$greencard[m3$race == "Pre_Black"], m3$greencard[m3$race ==    "Post_Black"], eps_std = 0.01) %>% summary() 
equiv.t.test(m3$greencard[m3$race == "Pre_Middle Eastern"], m3$greencard[m3$race == "Post_Middle Eastern"], eps_std = 0.24) %>% summary()   
equiv.t.test(m3$greencard[m3$race == "Pre_South Asian"], m3$greencard[m3$race == "Post_South Asian"], eps_std = 0.32) %>% summary() 

equiv.t.test(m3$greencard[m3$religion == "Pre_Christian"], m3$greencard[m3$religion == "Two-wave_Christian"], eps_std = 0.37) %>% summary() 
equiv.t.test(m3$greencard[m3$religion == "Pre_Jewish"], m3$greencard[m3$religion == "Two-wave_Jewish"], eps_std = 0.41) %>% summary() 
equiv.t.test(m3$greencard[m3$religion == "Pre_Muslim"], m3$greencard[m3$religion == "Two-wave_Muslim"], eps_std = 0.21) %>% summary() 

equiv.t.test(m3$greencard[m3$race == "Pre_White"], m3$greencard[m3$race == "Two-wave_White"], eps_std = 0.32) %>% summary() 
equiv.t.test(m3$greencard[m3$race == "Pre_Black"], m3$greencard[m3$race == "Two-wave_Black"], eps_std = 0.44) %>% summary() 
equiv.t.test(m3$greencard[m3$race == "Pre_Middle Eastern"], m3$greencard[m3$race == "Two-wave_Middle Eastern"], eps_std = 0.22) %>% summary()    
equiv.t.test(m3$greencard[m3$race == "Pre_South Asian"], m3$greencard[m3$race == "Two-wave_South Asian"], eps_std = 0.34) %>% summary() 

equiv.t.test(m3$greencard[m3$religion == "Post_Christian"], m3$greencard[m3$religion == "Two-wave_Christian"], eps_std = 0.41) %>% summary() 
equiv.t.test(m3$greencard[m3$religion == "Post_Jewish"], m3$greencard[m3$religion == "Two-wave_Jewish"], eps_std = 0.23) %>% summary() 
equiv.t.test(m3$greencard[m3$religion == "Post_Muslim"], m3$greencard[m3$religion == "Two-wave_Muslim"]  , eps_std = 0.37) %>% summary() 

equiv.t.test(m3$greencard[m3$race == "Post_White"], m3$greencard[m3$race == "Two-wave_White"], eps_std = 0.3) %>% summary() 
equiv.t.test(m3$greencard[m3$race == "Post_Black"], m3$greencard[m3$race == "Two-wave_Black"], eps_std = 0.45) %>% summary() 
equiv.t.test(m3$greencard[m3$race == "Post_Middle Eastern"], m3$greencard[m3$race == "Two-wave_Middle Eastern"], eps_std = 0.28) %>% summary() 
equiv.t.test(m3$greencard[m3$race == "Post_South Asian"], m3$greencard[m3$race == "Two-wave_South Asian"], eps_std = 0.22) %>% summary() 

```
